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Abstract 

In previous papers [Phys. Rev. A 41, 4501 (1990), Phys. Rev. E 18, 3178 
(1993)], simple equilibrium expressions were obtained for nonlinear Burnett 
coefficients. A preliminary calculation of a 32 particle Lennard- Jones fluid 
was presented in the previous paper. Now, sufficient resources have become 
available to address the question of whether nonlinear Burnett coefficients are 
finite for soft spheres. The hard sphere case is known to have infinite non- 
linear Burnett coefficients (ie a nonanalytic constitutive relation) from mode 
coupling theory. This paper reports a molecular dynamics caclulation of the 
third order nonlinear Burnett coefficient of a Lennard-Jones fluid undergoing 
colour flow, which indicates that this term is diverges in the thermodynamic 
limit. 
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I. INTRODUCTION 



Ever since the Green-Kubo formalism for calculating the linear transport coefficients 
was developed, there has been interest in a corresponding theory for the nonlinear Burnett 
coefficients. The discovery of long-time tails in the velocity autocorrelation function by Alder 
and Wainwright |l] indicated that the hydrodynamic transport coefficients do not exist in 
two dimensions, but do exist in three dimensions. By applying mode-coupling theories, Ernst 
et al. showed that the relation between stress and strain rate should be P xy oc I7I In I7I 
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for hard disks and P xy = —777 + cyyyi for hard spheres, which are non-analytic constitutive 
relations. Similar results were obtained by Kawasaki and Gunton || for incompressible fluids 
(which is a particular case of a hard sphere fluid), although criticised later by Brey et. al. 
0]. It should be pointed out that the linear Burnett coefficients are known to be divergent 
PH, and in light of the linear coefficient results, it is generally assumed that the nonlinear 
coefficients should be divergent as well for soft particle systems. Brey et. al. |4| claim to 
show the divergence of nonlinear coefficients in a followup paper (ref. 11 in that paper), 
yet this paper never appeared in the literature. Therefore, there is considerable interest in 
a molecular dynamics simulation of a soft particle system to see if the hard sphere results 
generalise. 

In a paper by Evans and Lynden-Bell , equilibrium fluctuation expressions for inverse 
Burnett coefficients were derived for the colour conductivity problem. The coefficients, 
E>i, give a Taylor series representation of a nonlinear transport coefficient L, in terms of 
the thermodynamic force F. Thus if a thermodynamic flux J is written in terms of the 
coefficient's defining constitutive relation as (J) = L(F)F, then the Burnett coefficients are 
related by L(F) = Bq + B\F + B2F 2 + ■ ■ -. In order to derive closed form expressions for the 
Burnett coefficients, it was found necessary to work in the Norton ensemble, in which the flux 
J, rather than the thermodynamic force F was the independent variable. The constitutive 
relation in this case is (F) = £( J) J = B + B\ J + ■ ■ -. In the thermodynamic limit, we 
may write C(J) = L _1 (J), and so the non-linear Burnett coefficients can be computed by 
inverting the series. 

Evans and Lynden-Bell applied constant current dynamics to a canonical ensemble 
with the currents distributed about an average current J . This allowed the derivation of a 
transient time correlation function for the non-equilibrium phase average (F). It was then a 
simple matter to compute the derivatives of (F) with respect to the average current Jo, as the 
constant current propagator commutes with the derivative operator. However, this method 
appeared to be limited to colour currents, for which an appropriate canonical distribution 
could be found. In a previous paper jSj we show that this method can be applied to the 
situation of an arbitrary thermodynamic flux. Later, || we showed that this transient time 
correlation expression can be expressed in terms of an average over an equilibrium simulation, 
reducing the calculation required by two orders of magnitude. At the time, computational 
resources were not sufficient to establish whether this expression is finite in the limit as 
t — > 00, or in the thermodynamic limit. In this paper, we present computational results of 
colour conductivity in a Lennard- Jones system, harnessing 4 supercomputers simultaneously 
over a period of 18 months, that show distinct evidence that B2 = 00. 

In order to avoid confusion, it should be noted that the term "colour diffusion" is some- 
times used in the sense of the diffusion of colour labels attached to otherwise colour blind 
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particles in the complete absence of applied external fields [|Kj . In this approach if the colour 



label attached to a particle is ignored, the system remains at equilibrium. This is manifestly 
a linear process. In the model we consider all the particles interact with an external colour 
sensitive external field and this allows the possibility of a nonlinear response. It might also 
be pointed out the the colour field we consider here is independent of both position and time 
so that the linear Burnett coefficients do not play a role. 



II. THE SIMULATION 

The simulation was performed using the colour conductivity model described in Evans 
and Lynden-Bell 0. The intermolecular potential was taken to be the Lennard- Jones po- 
tential, which has an attractive component due to van der Waals interaction, and a repulsive 
hard core that goes as r~ 12 : 

V{r) = 

In what follows, every quantity will be given in reduced units, in which e = a = m = 1. 
This model has been well studied, and can be related physically to some noble gases like 
argon. 

The system was simulated at 3 different system sizes (32, 108 and 256 particles) using 
a periodic boundary condition to minimise boundary effects. The state point chosen had a 
temperature of 1.08 and density of 0.85. Considerable information was already known about 




this system at that state point [11 



The equations of motion are just that of the Nose- Hoover thermostat, with an additional 
flux statting term. This generates a canonical ensemble: 

Pi 

= — , 

m 

pi = Fj + e^xA - api, 

a = ^r (T " To) ' 

X = ^(J-J(t = 0)), (1) 
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where Fj are the intermolecular forces, e« = ±1 are the colour charges, T = ^Nk B ' = ^ 
and J = J2i ^Vn7 * s ^ ne colour current. 

The feedback parameter Q\ was chosen equal to 4.74 for the 108, 256 and one of the 32 
particle runs. Because Q\ should be an extensive quantity, the 32 particle run was repeated 
at Q\ = 32 x 4.74/108 = 1.4. The Nose-Hoover thermostat parameter Q a was chosen to be 
0.3L/V. The values of these parameters were chosen to give optimal convergence of the linear 
response function. There is no real reason for them to be optimal for non-linear response 
functions. 

When the flux is fixed in this manner, the ensemble is termed a Norton ensemble. When 
the thermodynamic force is fixed, then it is termed a Thevenin ensemble by analogy with 
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electrical circuits We have recently given a statistical mechanical proof of the macro- 



scopic equivalence of the Norton and Thevenin representations of a nonequilibrium system 
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Recall that Transient Time Correlation Functions for evaluating the inverse cubic Burnett 
coefficient £>2 is given in Ref. J7|: 

62 = <tj^ / °°( A (^)A(0)(AJ 2 - <AJ 2 »)cfe. (2) 

where A(s) is the additional phase variable (defined in eq. (|I|) corresponding to a colour 
force of a system at time s along a trajectory and J is the colour current at the origin 
of that trajectory. As the system is at equilibrium (in the canonical ensemble), after a 
correlation time has passed, the system's configuration is effectively randomised, and may 
be used as a new trajectory origin. The correlations between different successive states 
of the equilibrium simulation can be easily seen by examining something like the velocity 
autocorrelation function (see Fig 7.1 of [|HJ for examples). The correlation time for this 
system is about 1. 



III. RESULTS 

Because the relevant quantity is an ensemble average, a very effective parallelisation 
strategy is to run a separate copy of the system on each processor, compute the TTCF on 
each processor, then average over the entire set of processors, weighting for the number of 
timesteps executed on each processor. Further computational details of this experiment have 
been reported in [|14]]. Whilst the results of this experiment would appear meagre compared 
with the computational resources used to compute it, it should be pointed out that this 
computation was conducted at the lowest priority on these machines, using idle CPU cycles. 

Having a set of approximations also allows one to calculate the standard error of the 
TTCF. These are shown as error bars in Figures 

The TTCFs and their integrals are shown in figures [|^. There is a considerable system 
size dependence, indicating that the nonlinear Burnett coefficients diverge in the thermo- 
dynamic limit, although the individual TTCFs remain finite. It can be shown, using the 
lemma proved in the appendix of f9j, that the inverse nonlinear Burnett coefficients given 
by equation (|2[) should be intensive. As well as this, the 32 particle simulation shows strong 
evidence of a long time tail (Fig ^| and |3|) when Q\ is increased (softening the current- 
statting), leading to a divergence in the integrals as t — > 00. For comparison, the transient 
time correlation function for the linear coefficient is shown in Fig. [I], showing convergence 
within t = 5. 



IV. CONCLUSION 

This work presents strong numerical evidence in favour of infinite nonlinear Burnett 
coefficients for soft spheres as is the case for hard spheres. However, the Taylor series 
expansion of the constitutive relation presented in || can also be derived for J 7^ 0, which 
if the hard sphere model is anything to go by, should be finite. These can be used to 
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compute the constitutive relation into the nonlinear region. However, it will probably be at 
least another decade before these calculations become practical. 
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FIG. 3. Integral of TTCF for the 32 particle system with Q x = 4.74 at 1.32 x 10 11 timesteps 




FIG. 4. Transient Time Correlation Function for the 32 particle system with Qx = 1.4 at 
2.2 x 10 11 timesteps 
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